写在前面

2023年01月28日,更新ggClusterNet文档。

R包安装

导入R包

#--导入所需R包#-------
library(phyloseq)
library(igraph)
library(network)
library(sna)
library(tidyverse)
library(ggClusterNet)

输入方式

输入方式一

直接输入phyloseq格式的数据。

#-----导入数据#-------
data(ps)

输入方式二

可以从https://github.com/taowenmicro/R-_function下载数据,构造phylsoeq文件。自己的数据也按照网址示例数据进行准备。虽然phylsoeq对象不易用常规手段处理,但是组学数据由于数据量比较大,数据注释内容多样化,所以往往使用诸如phyloseq这类对象进行处理,并简化数据处理过程。ggClusterNet同样使用了phyloseq对象作为微生物网络的分析。

phyloseq对象构建过程如下,网络分析主要用到otu表格,后续pipeline流程可能用到分组文件metadata,如果按照分类水平山色或者区分模块则需要taxonomy。这几个部分并不是都必须加入phyloseq对象中,可以用到那个加那个。

或者直接从网站读取,但是由于github在国外,所以容易失败

使用ggCLusterNet进行网络分析的过程

corMicro函数用于计算网络相关

按照丰度过滤微生物表格,并却计算相关矩阵,按照指定的阈值挑选矩阵中展示的数值。调用了psych包中的corr.test函数,使用三种相关方法。 N参数提取丰度最高的150个OTU; method.scale参数确定微生物组数据的标准化方式,这里我们选用TMM方法标准化微生物数据。

#-提取丰度最高的指定数量的otu进行构建网络


#----------计算相关#----
result = corMicro(ps = ps,
                   N = 150,
                   method.scale = "TMM",
                   r.threshold=0.8,
                   p.threshold=0.05,
                   method = "spearman"
                   
                   
                   )

#--提取相关矩阵
cor = result[[1]]
cor[1:6,1:6]
##         ASV_1 ASV_28 ASV_125 ASV_135 ASV_8 ASV_106
## ASV_1       1      0       0       0     0       0
## ASV_28      0      1       0       0     0       0
## ASV_125     0      0       1       0     0       0
## ASV_135     0      0       0       1     0       0
## ASV_8       0      0       0       0     1       0
## ASV_106     0      0       0       0     0       1
#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]

#-导出otu表格
otu_table = ps_net %>% 
  vegan_otu() %>%
  t() %>%
  as.data.frame()

制作分组,我们模拟五个分组

这是网络布局的基础,无论是什么聚类布局,都需要制作一个分组文件,这个文件有两列,一列是节点,一列是分组信息,这个分组信息名称为:group。这个文件信息就是用于对节点进行分组,然后按照分组对节点归类,使用包中可视化函数计算节点位置。

注意分组文件的格式,分为两列,第一列是网络中包含的OTU的名字,第二列是分组信息,同样的分组标记同样的字符。

#--人工构造分组信息:将网络中全部OTU分为五个部分,等分
netClu = data.frame(ID = row.names(otu_table),group =rep(1:5,length(row.names(otu_table)))[1:length(row.names(otu_table))] )
netClu$group = as.factor(netClu$group)
head(netClu)
##        ID group
## 1   ASV_1     1
## 2  ASV_28     2
## 3 ASV_125     3
## 4 ASV_135     4
## 5   ASV_8     5
## 6 ASV_106     1

PolygonClusterG 根据分组,计算布局位置坐标

不同的模块按照分组聚集成不同的圆,并且圆形的大小一样。如果一个分组只有一个点,则这个点坐落在圆心位置。

#--------计算布局#---------
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu  )
node = result2[[1]]
head(node)
##                    X1   X2 elements
## ASV_1    0.000000e+00 10.5    ASV_1
## ASV_14   2.598076e+00  9.0   ASV_14
## ASV_72   2.598076e+00  6.0   ASV_72
## ASV_109  3.673819e-16  4.5  ASV_109
## ASV_71  -2.598076e+00  6.0   ASV_71
## ASV_88  -2.598076e+00  9.0   ASV_88

nodeadd 节点注释的:用otu表格和分组文件进行注释

nodeadd函数只是提供了简单的用注释函数,用户可以自己在node的表格后面添加各种注释信息。

tax_table = ps_net %>%
  vegan_tax() %>%
  as.data.frame()
#---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
nodes[1:6,1:6]
##                X1         X2 elements  Kingdom         Phylum
## ASV_1    0.000000 10.5000000    ASV_1 Bacteria Actinobacteria
## ASV_10   3.792382 -8.9964485   ASV_10 Bacteria Proteobacteria
## ASV_100  7.269286 -5.1349547  ASV_100 Bacteria Proteobacteria
## ASV_101  5.911236  5.0628075  ASV_101 Bacteria Proteobacteria
## ASV_102  4.278276  1.3951201  ASV_102 Bacteria Proteobacteria
## ASV_103 -8.359017 -0.4411929  ASV_103 Bacteria Proteobacteria
##                       Class
## ASV_1        Actinobacteria
## ASV_10  Alphaproteobacteria
## ASV_100  Betaproteobacteria
## ASV_101          Unassigned
## ASV_102 Gammaproteobacteria
## ASV_103 Alphaproteobacteria

计算边

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
head(edge)
##         X2         Y2   OTU_2   OTU_1     weight        X1         Y1 cor
## 1 2.598076  9.0000000  ASV_14 ASV_172 -0.8020673  4.278276  3.2492221   -
## 2 7.131446  5.3221711  ASV_28  ASV_20  0.8041345  7.755181 -0.6122717   +
## 3 4.533370  0.8221711 ASV_140  ASV_64  0.8484250  1.818041 -4.5620057   +
## 4 7.014193 -4.5620057 ASV_166  ASV_43 -0.8064018 -2.853170  8.4270510   -
## 5 7.014193 -4.5620057 ASV_166 ASV_172 -0.8016529  4.278276  3.2492221   -
## 6 1.818041 -7.5620057  ASV_29  ASV_43 -0.8415076 -2.853170  8.4270510   -

出图

p <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank()) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p

# ggsave("cs1.pdf",p,width = 8,height = 6)

模拟不同的分组–可视化

模拟不同分组效果展示:3个分组

这是网络布局的基础,无论是什么聚类布局,都需要制作一个分组文件,这个文件有两列,一列是节点,一列是分组信息,这个分组信息名称必须设定为:group。

netClu = data.frame(ID = row.names(tax_table),group =rep(1,length(row.names(tax_table)))[1:length(row.names(tax_table))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs2.pdf",pnet,width = 7,height = 5.5)

模拟不同的分组查看效果:8个分组

netClu = data.frame(ID = row.names(cor),group =rep(1:8,length(row.names(cor)))[1:length(row.names(cor))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs3.pdf",pnet,width = 6,height = 5)

微生物分类分组可视化

#----------计算相关#----
result = corMicro (ps = ps,
                   N = 200,
                   method.scale = "TMM",
                   r.threshold=0.8,
                   p.threshold=0.05,
                   method = "spearman"
                   )
#--提取相关矩阵
cor = result[[1]]
# head(cor)
#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]
#-导出otu表格
otu_table = ps_net %>% 
  vegan_otu() %>%
  t() %>%
  as.data.frame()
tax = ps_net %>% vegan_tax() %>%
  as.data.frame()
tax$filed = tax$Phylum
group2 <- data.frame(ID = row.names(tax),group = tax$Phylum)
group2$group  =as.factor(group2$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs3.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化1-圆环半径调整PolygonRrClusterG

结果发现这些高丰度OTU大部分属于放线菌门和变形菌门,其他比较少。所以下面我们按照OTU数量的多少,对每个模块的大小进行重新调整。

result2 = PolygonRrClusterG(cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs4.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化2-model_filled_circle

用实心点作为每个模块的布局方式

set.seed(12)
#-实心圆2
result2 = model_filled_circle(cor = cor,
                              culxy =TRUE,
                              da = NULL,# 数据框,包含x,和y列
                              nodeGroup = group2,
                              mi.size = 1,# 最小圆圈的半径,越大半径越大
                              zoom = 0.3# 不同模块之间距离
                              )
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs5.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化3 model_maptree_group

用实心点作为每个模块布局方式2:model_maptree_group,智能布局不同分组之间的距离,在美学上特征更明显一点。

set.seed(12)
#-实心圆2
result2 = model_maptree_group(cor = cor,
                              nodeGroup = group2,
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs6.pdf",pnet,width = 6,height = 5)

按照网络模块化分组

模块布局算法 model_maptree_group

按照物种组成分类完成网络分析其实并不常用,更多的是按照模块分组,进行网络可视化。

#--modulGroup函数用于计算模块并整理成分组信息
netClu  = modulGroup( cor = cor,cut = NULL,method = "cluster_fast_greedy" )

result2 = model_maptree_group(cor = cor,
                              nodeGroup = group2,
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]

# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
# head(nodes)

nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs7.pdf",pnet,width = 6,height = 5)

模块布局算法 model_maptree2

使用升级的model_maptree2:不在可以将每个模块独立区分,而是将模块聚拢,并在整体布局上将离散的点同这些模块一同绘制到同心圆内。控制了图形的整体布局为圆形。

result2 = model_maptree2(cor = cor,
                              method = "cluster_fast_greedy"
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]

# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)
##                  X1           X2 elements  Kingdom         Phylum
## ASV_1   -13.8901865  -0.03852081    ASV_1 Bacteria Actinobacteria
## ASV_10   -0.4097123   3.33964715   ASV_10 Bacteria Proteobacteria
## ASV_100  -1.1227948 -12.61245151  ASV_100 Bacteria Proteobacteria
## ASV_101   5.7161404  20.40310410  ASV_101 Bacteria Proteobacteria
## ASV_102 -18.4899174 -10.64072413  ASV_102 Bacteria Proteobacteria
## ASV_103  16.2557654  10.33992920  ASV_103 Bacteria Proteobacteria
##                       Class           Order              Family           Genus
## ASV_1        Actinobacteria Actinomycetales Thermomonosporaceae      Unassigned
## ASV_10  Alphaproteobacteria     Rhizobiales        Rhizobiaceae       Rhizobium
## ASV_100  Betaproteobacteria Burkholderiales      Comamonadaceae  Hydrogenophaga
## ASV_101          Unassigned      Unassigned          Unassigned      Unassigned
## ASV_102 Gammaproteobacteria      Unassigned          Unassigned      Unassigned
## ASV_103 Alphaproteobacteria     Rhizobiales  Phyllobacteriaceae Phyllobacterium
##                              Species       mean
## ASV_1                     Unassigned 1473.61111
## ASV_10                    Unassigned  428.27778
## ASV_100    Hydrogenophaga_intermedia   61.00000
## ASV_101                   Unassigned   60.72222
## ASV_102                   Unassigned   64.00000
## ASV_103 Phyllobacterium_bourgognense   53.77778
nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

# ggsave("cs8.pdf",pnet,width = 6,height = 5)

上千OTU相关性计算测试model_igraph2布局

这个布局最近几年文章上使用非常多。

# cor_Big_micro2 增加了标准化方法和p值矫正方法
result = cor_Big_micro2(ps = ps,
                       N = 1000,
                       r.threshold=0.85,
                       p.threshold=0.05,
                       method = "pearson",
                       scale = FALSE
)

#--提取相关矩阵
cor = result[[1]]
dim(cor)
## [1] 1000 1000
# model_igraph2
result2 <- model_igraph2(cor = cor,
                         method = "cluster_fast_greedy",
                         seed = 12
)
node = result2[[1]]
dim(node)
## [1] 293   3
dat = result2[[2]]
head(dat)
##   orig_model      model   color      OTU        X1         X2
## 1         45 mini_model #C1C1C1 ASV_1591  7.835825  12.339063
## 2         15   model_15 #50A9AF  ASV_688 -1.523618  -9.873582
## 3         34   model_34 #FEF6B1    ASV_8  4.877237  13.710152
## 4          6    model_6 #F99655  ASV_174  2.811909 -10.288624
## 5          6    model_6 #F99655  ASV_716  4.314573 -10.310663
## 6         34   model_34 #FEF6B1  ASV_106  3.768050  13.839837
tem = data.frame(mod = dat$model,col = dat$color) %>%  
  dplyr::distinct( mod, .keep_all = TRUE)  
col = tem$col
names(col) = tem$mod

#---node节点注释#-----------
otu_table = as.data.frame(t(vegan_otu(ps)))
tax_table = as.data.frame(vegan_tax(ps))
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)
##                 X1         X2 elements  Kingdom         Phylum
## ASV_10   -2.259642  10.577060   ASV_10 Bacteria Proteobacteria
## ASV_100  -7.994668   8.125110  ASV_100 Bacteria Proteobacteria
## ASV_102   6.730077 -13.330438  ASV_102 Bacteria Proteobacteria
## ASV_104  -2.682913  12.612126  ASV_104 Bacteria Proteobacteria
## ASV_1043  3.463727  -9.298511 ASV_1043 Bacteria  Acidobacteria
## ASV_1048 12.307434  -7.316681 ASV_1048 Bacteria Actinobacteria
##                        Class           Order            Family          Genus
## ASV_10   Alphaproteobacteria     Rhizobiales      Rhizobiaceae      Rhizobium
## ASV_100   Betaproteobacteria Burkholderiales    Comamonadaceae Hydrogenophaga
## ASV_102  Gammaproteobacteria      Unassigned        Unassigned     Unassigned
## ASV_104   Betaproteobacteria Burkholderiales    Comamonadaceae    Polaromonas
## ASV_1043  Acidobacteria_Gp10      Unassigned        Unassigned           Gp10
## ASV_1048      Actinobacteria Actinomycetales Microbacteriaceae      Agromyces
##                            Species       mean
## ASV_10                  Unassigned 428.277778
## ASV_100  Hydrogenophaga_intermedia  61.000000
## ASV_102                 Unassigned  64.000000
## ASV_104                 Unassigned  60.444444
## ASV_1043                Unassigned   4.833333
## ASV_1048         Agromyces_indicus   4.722222
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
colnames(edge)[8] = "cor"
head(edge)
##          X2        Y2    OTU_2   OTU_1    weight        X1         Y1 cor
## 1  7.835825 12.339063 ASV_1591 ASV_710 0.9060991  7.958023  13.239291   +
## 2 -1.523618 -9.873582  ASV_688 ASV_596 0.9262611 -2.153803 -10.584342   +
## 3 -1.523618 -9.873582  ASV_688 ASV_780 0.8884699 -1.365654 -10.865238   +
## 4 -1.523618 -9.873582  ASV_688  ASV_20 0.8557562 -1.207825  -8.872789   +
## 5  4.877237 13.710152    ASV_8 ASV_106 0.8734179  3.768050  13.839837   +
## 6  4.877237 13.710152    ASV_8 ASV_334 0.9389280  4.640294  14.523675   +
tem2 = dat %>% 
  dplyr::select(OTU,model,color) %>%
  dplyr::right_join(edge,by =c("OTU" = "OTU_1" ) ) %>%
  dplyr::rename(OTU_1 = OTU,model1 = model,color1 = color)
head(tem2)
##      OTU_1     model1  color1        X2         Y2   OTU_2    weight        X1
## 1 ASV_1591 mini_model #C1C1C1  7.958023  13.239291 ASV_710 0.9060991  7.835825
## 2  ASV_688   model_15 #50A9AF -2.153803 -10.584342 ASV_596 0.9262611 -1.523618
## 3  ASV_688   model_15 #50A9AF -1.365654 -10.865238 ASV_780 0.8884699 -1.523618
## 4  ASV_688   model_15 #50A9AF -1.207825  -8.872789  ASV_20 0.8557562 -1.523618
## 5    ASV_8   model_34 #FEF6B1  3.768050  13.839837 ASV_106 0.8734179  4.877237
## 6    ASV_8   model_34 #FEF6B1  4.640294  14.523675 ASV_334 0.9389280  4.877237
##          Y1 cor
## 1 12.339063   +
## 2 -9.873582   +
## 3 -9.873582   +
## 4 -9.873582   +
## 5 13.710152   +
## 6 13.710152   +
tem3 = dat %>% 
  dplyr::select(OTU,model,color) %>%
  dplyr::right_join(edge,by =c("OTU" = "OTU_2" ) ) %>%
  dplyr::rename(OTU_2 = OTU,model2 = model,color2 = color)
head(tem3)
##      OTU_2     model2  color2        X2        Y2   OTU_1    weight        X1
## 1 ASV_1591 mini_model #C1C1C1  7.835825 12.339063 ASV_710 0.9060991  7.958023
## 2  ASV_688   model_15 #50A9AF -1.523618 -9.873582 ASV_596 0.9262611 -2.153803
## 3  ASV_688   model_15 #50A9AF -1.523618 -9.873582 ASV_780 0.8884699 -1.365654
## 4  ASV_688   model_15 #50A9AF -1.523618 -9.873582  ASV_20 0.8557562 -1.207825
## 5    ASV_8   model_34 #FEF6B1  4.877237 13.710152 ASV_106 0.8734179  3.768050
## 6    ASV_8   model_34 #FEF6B1  4.877237 13.710152 ASV_334 0.9389280  4.640294
##           Y1 cor
## 1  13.239291   +
## 2 -10.584342   +
## 3 -10.865238   +
## 4  -8.872789   +
## 5  13.839837   +
## 6  14.523675   +
tem4 = tem2 %>%inner_join(tem3)
head(tem4)
##      OTU_1     model1  color1        X2         Y2   OTU_2    weight        X1
## 1 ASV_1591 mini_model #C1C1C1  7.958023  13.239291 ASV_710 0.9060991  7.835825
## 2  ASV_688   model_15 #50A9AF -2.153803 -10.584342 ASV_596 0.9262611 -1.523618
## 3  ASV_688   model_15 #50A9AF -1.365654 -10.865238 ASV_780 0.8884699 -1.523618
## 4  ASV_688   model_15 #50A9AF -1.207825  -8.872789  ASV_20 0.8557562 -1.523618
## 5    ASV_8   model_34 #FEF6B1  3.768050  13.839837 ASV_106 0.8734179  4.877237
## 6    ASV_8   model_34 #FEF6B1  4.640294  14.523675 ASV_334 0.9389280  4.877237
##          Y1 cor     model2  color2
## 1 12.339063   + mini_model #C1C1C1
## 2 -9.873582   +   model_15 #50A9AF
## 3 -9.873582   +   model_15 #50A9AF
## 4 -9.873582   +   model_15 #50A9AF
## 5 13.710152   +   model_34 #FEF6B1
## 6 13.710152   +   model_34 #FEF6B1
edge2 = tem4 %>% mutate(color = ifelse(model1 == model2,as.character(model1),"across"),
                        manual = ifelse(model1 == model2,as.character(color1),"#C1C1C1")
                        )


col_edge = edge2 %>% dplyr::distinct(color, .keep_all = TRUE)  %>% 
  select(color,manual)
col0 = col_edge$manual
names(col0) = col_edge$color

library(ggnewscale)

p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = color),
                              data = edge2, size = 1) +
  scale_colour_manual(values = col0) 

# ggsave("./cs1.pdf",p1,width = 16,height = 14)
p2 = p1 +
   new_scale_color() +
  geom_point(aes(X1, X2,color =model), data = dat,size = 4) +
  scale_colour_manual(values = col) +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank()) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2

# ggsave("./cs2.pdf",p2,width = 12,height = 10)

网络属性和节点属性

网络性质计算

22年6月升级后版本包括了16项网络属性,包括周集中老师21年NCC文章中全部属性

dat = net_properties(igraph)
head(dat)
##                      value
## num.edges      122.0000000
## num.pos.edges  101.0000000
## num.neg.edges   21.0000000
## num.vertices    49.0000000
## connectance      0.1037415
## average.degree   4.9795918
# 升级后包含的网络属性更多
dat = net_properties.2(igraph,n.hub = T)
head(dat,n = 16)
##                                               value
## num.edges(L)                            122.0000000
## num.pos.edges                           101.0000000
## num.neg.edges                            21.0000000
## num.vertices(n)                          49.0000000
## Connectance(edge_density)                 0.1037415
## average.degree(Average K)                 4.9795918
## average.path.length                       2.1710403
## diameter                                  4.7234210
## edge.connectivity                         0.0000000
## mean.clustering.coefficient(Average.CC)   0.5690439
## no.clusters                               5.0000000
## centralization.degree                     0.1879252
## centralization.betweenness                0.1533413
## centralization.closeness                  1.1468429
## RM(relative.modularity)                   0.2573947
## the.number.of.keystone.nodes              1.0000000
dat = net_properties.3(igraph,n.hub = T)
head(dat,n = 16)
##                                               value
## num.edges(L)                            122.0000000
## num.pos.edges                           101.0000000
## num.neg.edges                            21.0000000
## num.vertices(n)                          49.0000000
## Connectance(edge_density)                 0.1037415
## average.degree(Average K)                 4.9795918
## average.path.length                       2.1710403
## diameter                                  4.7234210
## edge.connectivity                         0.0000000
## mean.clustering.coefficient(Average.CC)   0.5690439
## no.clusters                               5.0000000
## centralization.degree                     0.1879252
## centralization.betweenness                0.1533413
## centralization.closeness                  1.1468429
## RM(relative.modularity)                  -8.4867310
## the.number.of.keystone.nodes              1.0000000
# 增加了网络模块性(modularity.net )和随机网络模块性(modularity_random )
# dat = net_properties.4(igraph,n.hub = F)
# head(dat,n = 16)

节点性质计算

nodepro = node_properties(igraph)
head(nodepro)
##        igraph.degree igraph.closeness igraph.betweenness igraph.cen.degree
## ASV_28             3        0.2569444           2.333333                 3
## ASV_48             1        0.2032967           0.000000                 1
## ASV_34             5        0.3274336         140.216667                 5
## ASV_18             4        0.2846154          20.916667                 4
## ASV_20             3        0.2569444          36.000000                 3
## ASV_44            11        0.3627451           7.874009                11

Zipi基于模块对OTU进行分类

result = cor_Big_micro2(ps = ps,
                       N = 500,
                       r.threshold=0.6,
                       p.threshold=0.05,
                       # method = "pearson",
                       scale = FALSE
)

#--提取相关矩阵
cor = result[[1]]

result4 = nodeEdge(cor = cor)
#提取变文件
edge = result4[[1]]
#--提取节点文件
node = result4[[2]]
igraph  = igraph::graph_from_data_frame(edge, directed = FALSE, vertices = node)
res = ZiPiPlot(igraph = igraph,method = "cluster_fast_greedy")
p <- res[[1]]
p

# ggsave("./cs2.pdf",p,width = 8,height = 7)

扩展-关键OTU挑选

Hub节点是在网络中与其他节点连接较多的节点,Hub微生物就是与其他微生物联系较为紧密的微生物,可以称之为关键微生物(keystone)

hub = hub_score(igraph)$vector %>%
  sort(decreasing = TRUE) %>%
  head(5) %>%
  as.data.frame()

colnames(hub) = "hub_sca"

ggplot(hub) +
  geom_bar(aes(x = hub_sca,y = reorder(row.names(hub),hub_sca)),stat = "identity",fill = "#4DAF4A")

# ggsave("./cs2.pdf",width = 5,height = 4)

对应随机网络构建和网络参数比对

result = random_Net_compate(igraph = igraph, type = "gnm", step = 100, netName = layout)
p1 = result[[1]]
sum_net = result[[4]]
p1

head(sum_net)
##                          KO        Means SD
## num.edges      1.078000e+03 1.078000e+03  0
## num.pos.edges  7.570000e+02 0.000000e+00  0
## num.neg.edges  3.210000e+02 0.000000e+00  0
## num.vertices   3.350000e+02 3.350000e+02  0
## connectance    1.926892e-02 1.926892e-02  0
## average.degree 6.435821e+00 6.435821e+00  0
# ggsave("./cs3.pdf",p1,width = 5,height = 4)

微生物网络流程

微生物组小网络:model_Gephi.2

使用network函数运行微生物网络全套分析:

  • 使用OTU数量建议少于250个,如果OTU数量为250个,同时计算zipi,整个运算过程为3-5min。
data("ps16s")
path = "./result_micro_200/"
dir.create(path)

result = network(ps = ps16s,
                 N = 200,
                layout_net = "model_Gephi.2",
                 r.threshold=0.6,
                 p.threshold=0.05,
                 label = FALSE,
                 path = path,
                 zipi = TRUE)
## [1] "Group1"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] "Group2"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] "Group3"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
# 多组网络绘制到一个面板
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]


plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 48,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 48,height = 16)


tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络:model_maptree2

大网络运算时间会比较长,这里我没有计算zipi,用时5min完成全部运行。 N=0,代表用全部的OTU进行计算。

  • 3000个OTU不计算zipi全套需要18min。
path = "./result_big_1000/"
dir.create(path)

result = network.2(ps = ps16s,
                 N = 1000,
                 big = TRUE,
                 maxnode = 5,
                 select_layout = TRUE,
                  layout_net = "model_maptree2",
                 r.threshold=0.4,
                 p.threshold=0.01,
                 label = FALSE,
                 path = path,
                 zipi = FALSE)
## [1] "Group1"
## [1] "cor matrix culculating over"
## [1] "Group2"
## [1] "cor matrix culculating over"
## [1] "Group3"
## [1] "cor matrix culculating over"
# 多组网络绘制到一个面板
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10*num,height = 10,limitsize = FALSE)

tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络:model_igraph

path = "./result_1000_igraph/"
dir.create(path)
# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map
result = network.2(ps = ps16s,
                 N = 1000,
                 big = TRUE,
                 maxnode = 5,
                 select_layout = TRUE,
                  layout_net = "model_igraph",
                 r.threshold=0.4,
                 p.threshold=0.01,
                 label = FALSE,
                 path = path,
                 zipi = FALSE)
## [1] "Group1"
## [1] "cor matrix culculating over"
## [1] "Group2"
## [1] "cor matrix culculating over"
## [1] "Group3"
## [1] "cor matrix culculating over"
# 多组网络绘制到一个面板
p = result[[1]]
# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 16*num,height = 16,limitsize = FALSE)

tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络:model_igraph2-network.i

model_igraph2在model_igraph的基础上去除了离散点,这对于重复数量少的研究比较有利,可以更加清楚的展示小样本网络。

network.i函数专门为model_igraph2算法而写,可以额外输出模块上色的网络。

path = "./result_1000_igraph2/"
dir.create(path)

# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map

result = network.i(ps =  ps16s,
                   N = 1000,
                   r.threshold=0.9,
                   big = T,
                   select_layout = T,
                   method = "pearson",
                   scale = FALSE,
                   layout_net = "model_igraph2",
                   p.threshold=0.05,
                   label = FALSE,
                   path =path ,
                   zipi = F,
                   order = NULL
                   )
## [1] "Group1"
## [1] "cor matrix culculating over"
## [1] "Group2"
## [1] "cor matrix culculating over"
## [1] "Group3"
## [1] "cor matrix culculating over"
p1 = result[[1]]
p1

dat = result[[2]]
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(dat,tablename)
p = result[[5]]


plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p1,width = 10,height = 3,limitsize = FALSE)

plotname1 = paste(path,"/network_all2.pdf",sep = "")
ggsave(plotname1, p,width = 25,height = 8,limitsize = FALSE)

跨域网络流程

细菌,真菌,环境因子三者跨域网络

#--细菌和环境因子的网络#--------
Envnetplot<- paste("./16s_ITS_Env_network",sep = "")
dir.create(Envnetplot)

ps16s =ps16s%>% ggClusterNet::scale_micro()
psITS = psITS%>% ggClusterNet::scale_micro()


#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS= psITS,
                         NITS = 200,
                         N16s = 200)

map =  phyloseq::sample_data(ps.merge)
head(map)
##                ID  Group
## sample1   sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map

#--环境因子导入
data1 = env
envRDA  = data.frame(row.names = env$ID,env[,-1])

envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s

Gru = data.frame(ID = colnames(env)[-1],group = "env" )
head(Gru)
##      ID group
## 1    pH   env
## 2   SOC   env
## 3    TN   env
## 4 NH4.N   env
## 5 NO3.N   env
## 6    AP   env
library(sna)
library(ggClusterNet)
library(igraph)

result <- ggClusterNet::corBionetwork(ps = ps.merge,
                        N = 0,
                        r.threshold = 0.4, # 相关阈值
                        p.threshold = 0.05,
                        big = T,
                        group = "Group",
                        env = data1, # 环境指标表格
                        envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 18,
                        label = TRUE,
                        height = 10
)
## [1] "one"
##  num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
##   ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...
## [1] "1"
## [1] "2"
## [1] "3"
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 15,height = 12,dpi = 72)
# plotname1 = paste(Envnetplot,"/network_all.png",sep = "")
# ggsave(plotname1, p,width = 10,height = 8,dpi = 72)
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 15,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

细菌-环境因子跨域网络

#--细菌-环境因子网络#-----
Envnetplot<- paste("./16s_Env_network",sep = "")
dir.create(Envnetplot)
# ps16s = ps0 %>% ggClusterNet::scale_micro()
psITS = NULL

#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = NULL,
                                       N16s = 500)
ps.merge
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 500 taxa and 18 samples ]
## sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 500 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)
head(map)
##                ID  Group
## sample1   sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
data(env)
data1 = env
envRDA  = data.frame(row.names = env$ID,env[,-1])

envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s

Gru = data.frame(ID = colnames(env)[-1],group = "env" )
head(Gru)
##      ID group
## 1    pH   env
## 2   SOC   env
## 3    TN   env
## 4 NH4.N   env
## 5 NO3.N   env
## 6    AP   env
# library(sna)
# library(ggClusterNet)
# library(igraph)

result <- ggClusterNet::corBionetwork(ps = ps.merge,
                                      N = 0,
                                      r.threshold = 0.4, # 相关阈值
                                      p.threshold = 0.05,
                                      big = T,
                                      group = "Group",
                                      env = data1, # 环境指标表格
                                      envGroup = Gru,# 环境因子分组文件表格
                                      # layout = "fruchtermanreingold",
                                      path = Envnetplot,# 结果文件存储路径
                                      fill = "Phylum", # 出图点填充颜色用什么值
                                      size = "igraph.degree", # 出图点大小用什么数据
                                      scale = TRUE, # 是否要进行相对丰度标准化
                                      bio = TRUE, # 是否做二分网络
                                      zipi = F, # 是否计算ZIPI
                                      step = 100, # 随机网络抽样的次数
                                      width = 18,
                                      label = TRUE,
                                      height = 10
)
## [1] "one"
##  num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
##   ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...
## [1] "1"
## [1] "2"
## [1] "3"
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 13,height = 12,dpi = 72)

plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 13,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

细菌和真菌跨域网络

#---细菌和真菌OTU网络-域网络-二分网络#-------
# 仅仅关注细菌和真菌之间的相关,不关注细菌内部和真菌内部相关

Envnetplot<- paste("./16S_ITS_network",sep = "")
dir.create(Envnetplot)


data(psITS)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = psITS,
                                       N16s = 500,
                                       NITS = 500
)

ps.merge
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1000 taxa and 18 samples ]
## sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 1000 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)

head(map)
##                ID  Group
## sample1   sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map


data =  data.frame(ID = c("fun_ASV_205","fun_ASV_316","fun_ASV_118"),
                   c("Coremode","Coremode","Coremode"))

# source("F:/Shared_Folder/Function_local/R_function/my_R_packages/ggClusterNet/R/corBionetwork2.R")


result <- corBionetwork(ps = ps.merge,
                        N = 0,
                        lab = data,
                        r.threshold = 0.8, # 相关阈值
                        p.threshold = 0.05,
                        group = "Group",
                        # env = data1, # 环境指标表格
                        # envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 12,
                        label = TRUE,
                        height = 10,
                        big = TRUE,
                        select_layout = TRUE,
                        layout_net = "model_maptree2",
                        clu_method = "cluster_fast_greedy"
                        
                        
)
## [1] "one"
## [1] "1"
## [1] "2"
## [1] "3"
tem <- model_maptree(cor =result[[5]],
                     method = "cluster_fast_greedy",
                     seed = 12
)
node_model = tem[[2]]
head(node_model)
##            ID group degree
## 1 fun_ASV_246    40     35
## 2 fun_ASV_770    40     33
## 3 fun_ASV_939    40     30
## 4 fun_ASV_233    41     22
## 5  fun_ASV_76    43     22
## 6 fun_ASV_268    19     19
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]

plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 20,height = 19)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)
tablename <- paste(Envnetplot,"/node_model_imformation",".csv",sep = "")
write.csv(node_model,tablename)

tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)

细菌真菌任意分类水平跨域网络

#--细菌真菌任意分类水平跨域网络#----------
library(tidyverse)
# res1path = "result_and_plot/Micro_and_other_index_16s/"
Envnetplot<- paste("./16S_ITS_network_Genus",sep = "")
dir.create(Envnetplot)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- merge16S_ITS(ps16s = ps16s,
                         psITS = psITS,
                         N16s = 500,
                         NITS = 500
)
ps.merge
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1000 taxa and 18 samples ]
## sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
## tax_table()   Taxonomy Table:    [ 1000 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)
head(map)
##                ID  Group
## sample1   sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
tem.0 = ps.merge %>% tax_glom_wt(ranks = "Genus")
tax = tem.0 %>% vegan_tax() %>%
  as.data.frame()
head(tax)
##                    filed  Kingdom         Phylum               Class
## Acaulium             fun    Fungi     Ascomycota     Sordariomycetes
## Acidocella           bac Bacteria Proteobacteria Alphaproteobacteria
## Acrobeloides         fun  Metazoa       Nematoda         Chromadorea
## Actinomadura         bac Bacteria Actinobacteria      Actinobacteria
## Altererythrobacter   bac Bacteria Proteobacteria Alphaproteobacteria
## Alternaria           fun    Fungi     Ascomycota     Dothideomycetes
##                               Order              Family              Genus
## Acaulium               Microascales       Microascaceae           Acaulium
## Acidocella         Rhodospirillales    Acetobacteraceae         Acidocella
## Acrobeloides             Rhabditida        Cephalobidae       Acrobeloides
## Actinomadura        Actinomycetales Thermomonosporaceae       Actinomadura
## Altererythrobacter Sphingomonadales  Erythrobacteraceae Altererythrobacter
## Alternaria             Pleosporales       Pleosporaceae         Alternaria
data =  data.frame(ID = c("Acaulium","Acidocella","Acrobeloides"),
                   c("Acaulium","Acidocella","Acrobeloides"))

result <- corBionetwork(ps = tem.0,
                        N = 0,
                        lab = data,
                        r.threshold = 0.6, # 相关阈值
                        p.threshold = 0.05,
                        group = "Group",
                        # env = data1, # 环境指标表格
                        # envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 12,
                        label = TRUE,
                        height = 10,
                        big = TRUE,
                        select_layout = TRUE,
                        layout_net = "model_maptree2",
                        clu_method = "cluster_fast_greedy"
                        
                        
)
## [1] "one"
## [1] "1"
## [1] "2"
## [1] "3"
tem <- model_maptree(cor =result[[5]],
                     method = "cluster_fast_greedy",
                     seed = 12
)
node_model = tem[[2]]
head(node_model)
##               ID group degree
## 1     Microascus    12     13
## 2 Scopulariopsis    13     12
## 3      Arcopilus    10     10
## 4   Cladosporium    13      9
## 5       Conocybe    13      8
## 6    Hydrogonium     6      8
library(WGCNA)
otu = tem.0 %>% vegan_otu() %>%
  as.data.frame()
node_model = node_model[match(colnames(otu),node_model$ID),]


MEList = moduleEigengenes(otu, colors = node_model$group)
MEs = MEList$eigengenes

nGenes = ncol(otu)
nSamples = nrow(otu)
moduleTraitCor = cor(MEs, envRDA, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

#sizeGrWindow(10,6)
pdf(file=paste(Envnetplot,"/","Module-env_relationships.pdf",sep = ""),width=10,height=6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")

dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(envRDA),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = greenWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))
dev.off()
## png 
##   2
p = result[[1]]
p

# 全部样本网络参数比对
data = result[[2]]


plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10,height = 10)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)


tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)

网络稳定性(抗干扰性)

模块相似度

#--网络稳定性:模块相似性------
data(ps)

library(tidyfst)

res = module.compare.m(
    ps = ps,
    Top = 500,
    degree = TRUE,
    zipi = FALSE,
    r.threshold= 0.8,
    p.threshold=0.05,
    method = "spearman",
    padj = F,
    n = 3)
## [1] 80  8
## [1] 159   8
## [1] 224   8
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
#不同分组使用一个圆圈展示,圆圈内一个点代表一个模块,相连接的模块代表了相似的模块。
p = res[[1]]
p

 #--提取模块的OTU,分组等的对应信息
dat = res[[2]]
head(dat)
##        ID     group Group
## 1 ASV_256 KOmodel_1    KO
## 2 ASV_142 KOmodel_1    KO
## 3 ASV_362 KOmodel_1    KO
## 4  ASV_89 KOmodel_1    KO
## 5 ASV_419 KOmodel_1    KO
## 6 ASV_352 KOmodel_1    KO
#模块相似度结果表格
dat2 = res[[3]]
head(dat2)
##                          module1    module2 Both P1A2 P2A1 A1A2
## KOmodel_41-OEmodel_2  KOmodel_41  OEmodel_2    1    2    6  500
## KOmodel_62-OEmodel_1  KOmodel_62  OEmodel_1    1    2    6  500
## KOmodel_34-OEmodel_1  KOmodel_34  OEmodel_1    1    2    6  500
## KOmodel_22-OEmodel_4  KOmodel_22  OEmodel_4    1    3    5  500
## KOmodel_26-OEmodel_8  KOmodel_26  OEmodel_8    1    3    4  501
## KOmodel_17-OEmodel_10 KOmodel_17 OEmodel_10    1    4    4  500
##                                    p_raw              p_adj
## KOmodel_41-OEmodel_2  0.0407716775257314 0.0407716775257314
## KOmodel_62-OEmodel_1  0.0407716775257314 0.0407716775257314
## KOmodel_34-OEmodel_1  0.0407716775257314 0.0407716775257314
## KOmodel_22-OEmodel_4  0.0464588019672787 0.0464588019672787
## KOmodel_26-OEmodel_8  0.0388304723830171 0.0388304723830171
## KOmodel_17-OEmodel_10 0.0483470023594227 0.0483470023594227

网络鲁棒性(随机去除节点)

这里通过随机去除部分OTU,计算网络鲁棒性,代表网络抗干扰的能力。 按照0.05的步长,每次去除5%的文生物,重新计算鲁棒性,知道最终全部去除。 如果有分组列Group,则会按照分组进行鲁棒性计算,并且区分颜色绘制到一个面板中。 计算鲁棒性这里使用丰度加成权重和不加权两种方式,左边是不加权,后侧是加权的结果。 这里步长不可以修改,因为修改好像也没什么意思。

#--随即取出任意比例节点-网络鲁棒性#---------
res = Robustness.Random.removal(ps = ps,
                                Top = 500,
                                r.threshold= 0.8,
                                p.threshold=0.05,
                                method = "spearman"
                                
                                )


p = res[[1]]

p

#提取数据
dat = res[[2]]
head(dat)
##   Proportion.removed remain.mean  remain.sd   remain.se weighted Group
## 1               0.05   0.5408197 0.01310521 0.001310521 weighted    KO
## 2               0.10   0.5030601 0.01536428 0.001536428 weighted    KO
## 3               0.15   0.4674044 0.01755406 0.001755406 weighted    KO
## 4               0.20   0.4325683 0.01966370 0.001966370 weighted    KO
## 5               0.25   0.3944809 0.02034929 0.002034929 weighted    KO
## 6               0.30   0.3542623 0.02243074 0.002243074 weighted    KO

网络鲁棒性(去除关键节点)

#---去除关键节点-网络鲁棒性#------
res= Robustness.Targeted.removal(ps = ps,
                                 Top = 500,
                                 degree = TRUE,
                                 zipi = FALSE,
                                 r.threshold= 0.8,
                                 p.threshold=0.05,
                                 method = "spearman")
p = res[[1]]
p

#提取数据
dat = res[[2]]

负相关比例

#--计算负相关的比例#----
res = negative.correlation.ratio(ps = ps,
                          Top = 500,
                          degree = TRUE,
                          zipi = FALSE,
                          r.threshold= 0.8,
                          p.threshold=0.05,
                          method = "spearman")

p = res[[1]]
p

dat = res[[2]]
#-负相关比例数据
head(dat)
##   ID    ratio
## 1 KO 41.15139
## 2 OE 34.19023
## 3 WT 39.24051
# dir.create("./negative_correlation_ratio/")
# write.csv(dat,
#           paste("./negative_correlation_ratio/","_negative_ratio_network.csv",sep = ""))

组成稳定性

treat = ps %>% sample_data()
treat$pair = paste( "A",c(rep(1:6,3)),sep = "")
head(treat)
##     Group      Date    Site Sequencing  Platform     Species Batch
## KO1    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO2    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO3    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO4    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO5    KO  2017/7/4  Harbin        BGI HiSeq2500 Arabidopsis     1
## KO6    KO  2017/7/5  Harbin        BGI HiSeq2500 Arabidopsis     1
##     BarcodeSequence LinkerPrimerSequence      ReversePrimer  ID pair
## KO1      ACGCTCGACA  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO1   A1
## KO2      ATCAGACACG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO2   A2
## KO3      ATATCGCGAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO3   A3
## KO4      CACGAGACAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO4   A4
## KO5      CTCGCGTGTC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO5   A5
## KO6      TAGTATCAGC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO6   A6
sample_data(ps) = treat
#一般性的处理,没有时间梯度的,这里设定time为F,意味着每两个群落的组合都进行比较
res = community.stability( ps = ps,time = FALSE)
p = res[[1]]
p

dat = res[[2]]

# 如果是时间序列,会按照时间的单一流动性方向进行比较,自然就不是两两比对了。
res = community.stability( ps = ps,time = TRUE)
p = res[[1]]
p

网络抗毁性

网络抗毁性:使用了自然连通度的概念,用于反映网络稳定性.具体而言,通过在构建好的网络中移除里面的节点,并计算自然连通度。这样随着取出点的数量逐渐增多,就可以计算得到一连串的网络连通度。这个算法最先来自PNAS的文章:“Reduced microbial stability in the active layer is associated with carbon loss under alpine permafrost degradation”

library(tidyfst)
library("pulsar")
library(ggClusterNet)
library(phyloseq)
library(tidyverse)

res = natural.con.microp (
    ps = ps,
    Top = 200,
    r.threshold= 0.5,
    p.threshold=0.05,
    method = "spearman",
    norm = F,
    end = 150,# 小于网络包含的节点数量
    start = 0,
    con.method = "pulsar"
    )

p = res[[1]]
p

dat  = res[[2]]

网络模块化

展示网络模块

这里我们通过指定模块

library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)


resu  = module_display.2(
  pst = ps,
  r.threshold= 0.6,
  p.threshold=0.05,
  select.mod = c("model_1","model_2","model_3","model_4"),#选择指定模块可视化
  Top = 500,
  num = 5, # 模块包含OTU数量少于5个的不展示,
  leg.col = 3
)

# 全部模块输出展示
p1 = resu[[1]]
p1

#--提取率相关矩阵
dat0 = resu[[4]]

#--提取模块信息:这里区分好不同模块的OTU,以及每个模块包含的网络边的数量
dat = resu[[5]]
head(dat)
##        ID   group degree
## 1  ASV_66 model_3     46
## 2 ASV_294 model_3     40
## 3 ASV_326 model_3     34
## 4  ASV_60 model_2     31
## 5   ASV_2 model_3     31
## 6 ASV_163 model_3     30
#--提取网络的边和节点表格
dat2 = resu[[6]]
names(dat2)
## [1] "edge" "node"

指定的目标模块展示

p2 = resu[[2]]
p2

#按照模块着色,模块之间的边着色为灰色展示整个网络。
p2 = resu[[3]]
p2

#--注意这个报错信息,为展示的空间不够大,拉一拉展示框即可
# Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  :
#   Viewport has zero dimension(s)

基于模块的微生物组成分析

对微生物数据分好模块,使用模块分类的数据导入module_composition函数,即可作指定模块的物种组成。

#--这里我们指定三个模块,绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")
#--模块信息,共有三列,第一列是OTU,第二列是吗模块的分类,第三列是模块捏的边数量
mod1 = resu$mod.groups %>% filter(group %in% select.mod)
head(mod1)
##        ID   group degree
## 1  ASV_66 model_3     46
## 2 ASV_294 model_3     40
## 3 ASV_326 model_3     34
## 4  ASV_60 model_2     31
## 5   ASV_2 model_3     31
## 6 ASV_163 model_3     30
#-计算模块物种组成信息
pst = ps %>%
  filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
  scale_micro("rela") %>%
  # subset_samples.wt("Group" ,id3[m]) %>%
  filter_OTU_ps(500)

#-选择模块的微生物组成,通过参数j选择不同的分类水平。
res = module_composition(pst = pst,mod1 = mod1,j = "Family")
p3 = res[[1]]
p3

#按照属水平进行绘制
res = module_composition(pst = pst,mod1 = mod1,j = "Genus")
p3 = res[[1]]
p3

#---提取出图数据物种组成数据导出
ps.t = res[[3]]

otu = ps.t %>% vegan_otu() %>% t() %>%
  as.data.frame()
# head(otu)
tax = ps.t %>% vegan_tax() %>%
  as.data.frame()
# head(tax)
tab.res = cbind(otu,tax)
# head(tab.res)

# 没有标准化的出图数据
tab.res.2 = res[[4]]$bundance
# head(tab.res.2)
# 标准化到100%的出图数据
tab.res.3 = res[[4]]$relaabundance
# head(tab.res.3)

模块计算alpha多样性

这里将指定模块和ps对象输入函数module_alpha即可计算三种alpha多样性,这里调用了EasyStat包,使用非参数检验检测了alpha多样性的显著性。

注意:这里的ps对象需要原始count数据,不可用标准化后的数据。

#--这里我们指定三个模块,绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")

#--模块信息,两列,第一列是OTU,第二列是模块的分类
mod1 = resu$mod.groups %>% 
  dplyr::select(ID,group) %>%
  dplyr::filter(group %in% select.mod)
head(mod1)
##        ID   group
## 1  ASV_66 model_3
## 2 ASV_294 model_3
## 3 ASV_326 model_3
## 4  ASV_60 model_2
## 5   ASV_2 model_3
## 6 ASV_163 model_3
#-选择模块的多样性
result = module_alpha (
  ps = ps,
  mod1 = mod1)

p5 = result[[1]]
p5

#--alpha多样性数据
plotd = result[[4]]$alpha
# alpha多样性显著性表格
sigtab = result[[4]]$sigtab

基于单个样本计算网络属性

这里我们需要输入微生物组数据的ps对象,其次,需要输入这组微生物的相关矩阵;然后netproperties.sample函数会按照单个样本计算网络属性;这里一共有15中网络属性。

这里计算出来后便可以和其他指标进行相关性分析;

library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)
#-计算模块物种组成信息
pst = ps %>%
  filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
  scale_micro("rela") %>%
  # subset_samples.wt("Group" ,id3[m]) %>%
  filter_OTU_ps(500)

# 选择和网络属性做相关的指标
select.env = "env1"
#--内置数据无需导入

env = env1 %>% 
  as.data.frame() %>%
  rownames_to_column("id")
dat.f = netproperties.sample(pst = pst,cor = cor)
# head(dat.f)


dat.f$id = row.names(dat.f)
dat.f = dat.f %>% dplyr:: select(id,everything())


subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
tab = dat.f %>% left_join(subenv,by = "id")

mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"id"))
lab = mean(mtcars2[,select.env])
preoptab = tab
  
p0_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
    ggplot2::geom_point() +
    ggpubr::stat_cor(label.y=lab*1.1)+
    ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
    facet_wrap(~variable, scales="free_x") +
    geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
    theme_classic()

p0_1

计算模块的Zscore值

按照样本计算模块的ZS值,方便施用这一指标和其他相关指标进行相关性分析。

dat = ZS.net.module(
    pst = ps,
    Top = 500,
    r.threshold= 0.8,
    p.threshold=0.05,
    select.mod =  NULL
    )
# head(dat)


map =sample_data(ps)
map$id = row.names(map)
map = map[,c("id","Group")]
data = map %>%
    as.tibble() %>%
    inner_join(dat,by = "id") %>%
    dplyr::rename(group = Group)

colnames(env)[1] = "id"
subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
# head(data)
tab = data %>% left_join(subenv,by = "id")
modenv = tab
mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"group","id"))
mtcars2$variable
##   [1] model_1   model_1   model_1   model_1   model_1   model_1   model_1  
##   [8] model_1   model_1   model_1   model_1   model_1   model_1   model_1  
##  [15] model_1   model_1   model_1   model_1   model_3   model_3   model_3  
##  [22] model_3   model_3   model_3   model_3   model_3   model_3   model_3  
##  [29] model_3   model_3   model_3   model_3   model_3   model_3   model_3  
##  [36] model_3   model_2   model_2   model_2   model_2   model_2   model_2  
##  [43] model_2   model_2   model_2   model_2   model_2   model_2   model_2  
##  [50] model_2   model_2   model_2   model_2   model_2   model_4   model_4  
##  [57] model_4   model_4   model_4   model_4   model_4   model_4   model_4  
##  [64] model_4   model_4   model_4   model_4   model_4   model_4   model_4  
##  [71] model_4   model_4   model_7   model_7   model_7   model_7   model_7  
##  [78] model_7   model_7   model_7   model_7   model_7   model_7   model_7  
##  [85] model_7   model_7   model_7   model_7   model_7   model_7   model_31 
##  [92] model_31  model_31  model_31  model_31  model_31  model_31  model_31 
##  [99] model_31  model_31  model_31  model_31  model_31  model_31  model_31 
## [106] model_31  model_31  model_31  model_12  model_12  model_12  model_12 
## [113] model_12  model_12  model_12  model_12  model_12  model_12  model_12 
## [120] model_12  model_12  model_12  model_12  model_12  model_12  model_12 
## [127] model_5   model_5   model_5   model_5   model_5   model_5   model_5  
## [134] model_5   model_5   model_5   model_5   model_5   model_5   model_5  
## [141] model_5   model_5   model_5   model_5   model_23  model_23  model_23 
## [148] model_23  model_23  model_23  model_23  model_23  model_23  model_23 
## [155] model_23  model_23  model_23  model_23  model_23  model_23  model_23 
## [162] model_23  model_19  model_19  model_19  model_19  model_19  model_19 
## [169] model_19  model_19  model_19  model_19  model_19  model_19  model_19 
## [176] model_19  model_19  model_19  model_19  model_19  model_14  model_14 
## [183] model_14  model_14  model_14  model_14  model_14  model_14  model_14 
## [190] model_14  model_14  model_14  model_14  model_14  model_14  model_14 
## [197] model_14  model_14  model_6   model_6   model_6   model_6   model_6  
## [204] model_6   model_6   model_6   model_6   model_6   model_6   model_6  
## [211] model_6   model_6   model_6   model_6   model_6   model_6   model_15 
## [218] model_15  model_15  model_15  model_15  model_15  model_15  model_15 
## [225] model_15  model_15  model_15  model_15  model_15  model_15  model_15 
## [232] model_15  model_15  model_15  model_29  model_29  model_29  model_29 
## [239] model_29  model_29  model_29  model_29  model_29  model_29  model_29 
## [246] model_29  model_29  model_29  model_29  model_29  model_29  model_29 
## [253] model_13  model_13  model_13  model_13  model_13  model_13  model_13 
## [260] model_13  model_13  model_13  model_13  model_13  model_13  model_13 
## [267] model_13  model_13  model_13  model_13  model_8   model_8   model_8  
## [274] model_8   model_8   model_8   model_8   model_8   model_8   model_8  
## [281] model_8   model_8   model_8   model_8   model_8   model_8   model_8  
## [288] model_8   model_22  model_22  model_22  model_22  model_22  model_22 
## [295] model_22  model_22  model_22  model_22  model_22  model_22  model_22 
## [302] model_22  model_22  model_22  model_22  model_22  model_25  model_25 
## [309] model_25  model_25  model_25  model_25  model_25  model_25  model_25 
## [316] model_25  model_25  model_25  model_25  model_25  model_25  model_25 
## [323] model_25  model_25  model_30  model_30  model_30  model_30  model_30 
## [330] model_30  model_30  model_30  model_30  model_30  model_30  model_30 
## [337] model_30  model_30  model_30  model_30  model_30  model_30  model_26 
## [344] model_26  model_26  model_26  model_26  model_26  model_26  model_26 
## [351] model_26  model_26  model_26  model_26  model_26  model_26  model_26 
## [358] model_26  model_26  model_26  model_20  model_20  model_20  model_20 
## [365] model_20  model_20  model_20  model_20  model_20  model_20  model_20 
## [372] model_20  model_20  model_20  model_20  model_20  model_20  model_20 
## [379] model_9   model_9   model_9   model_9   model_9   model_9   model_9  
## [386] model_9   model_9   model_9   model_9   model_9   model_9   model_9  
## [393] model_9   model_9   model_9   model_9   model_16  model_16  model_16 
## [400] model_16  model_16  model_16  model_16  model_16  model_16  model_16 
## [407] model_16  model_16  model_16  model_16  model_16  model_16  model_16 
## [414] model_16  model_10  model_10  model_10  model_10  model_10  model_10 
## [421] model_10  model_10  model_10  model_10  model_10  model_10  model_10 
## [428] model_10  model_10  model_10  model_10  model_10  model_28  model_28 
## [435] model_28  model_28  model_28  model_28  model_28  model_28  model_28 
## [442] model_28  model_28  model_28  model_28  model_28  model_28  model_28 
## [449] model_28  model_28  model_33  model_33  model_33  model_33  model_33 
## [456] model_33  model_33  model_33  model_33  model_33  model_33  model_33 
## [463] model_33  model_33  model_33  model_33  model_33  model_33  model_18 
## [470] model_18  model_18  model_18  model_18  model_18  model_18  model_18 
## [477] model_18  model_18  model_18  model_18  model_18  model_18  model_18 
## [484] model_18  model_18  model_18  model_21  model_21  model_21  model_21 
## [491] model_21  model_21  model_21  model_21  model_21  model_21  model_21 
## [498] model_21  model_21  model_21  model_21  model_21  model_21  model_21 
## [505] model_17  model_17  model_17  model_17  model_17  model_17  model_17 
## [512] model_17  model_17  model_17  model_17  model_17  model_17  model_17 
## [519] model_17  model_17  model_17  model_17  model_11  model_11  model_11 
## [526] model_11  model_11  model_11  model_11  model_11  model_11  model_11 
## [533] model_11  model_11  model_11  model_11  model_11  model_11  model_11 
## [540] model_11  model_27  model_27  model_27  model_27  model_27  model_27 
## [547] model_27  model_27  model_27  model_27  model_27  model_27  model_27 
## [554] model_27  model_27  model_27  model_27  model_27  model_32  model_32 
## [561] model_32  model_32  model_32  model_32  model_32  model_32  model_32 
## [568] model_32  model_32  model_32  model_32  model_32  model_32  model_32 
## [575] model_32  model_32  model_24  model_24  model_24  model_24  model_24 
## [582] model_24  model_24  model_24  model_24  model_24  model_24  model_24 
## [589] model_24  model_24  model_24  model_24  model_24  model_24  mother_no
## [596] mother_no mother_no mother_no mother_no mother_no mother_no mother_no
## [603] mother_no mother_no mother_no mother_no mother_no mother_no mother_no
## [610] mother_no mother_no mother_no
## 34 Levels: model_1 model_3 model_2 model_4 model_7 model_31 ... mother_no
head(mtcars2)
##   env1 group  id variable      value
## 1 7.67    KO KO1  model_1  7.1549454
## 2 7.61    KO KO2  model_1 -0.3726327
## 3 7.15    KO KO3  model_1 10.6994573
## 4 6.89    KO KO4  model_1 10.3945282
## 5 7.39    KO KO5  model_1 35.8009267
## 6 7.01    KO KO6  model_1  9.5297211
lab = mean(mtcars2[,select.env])
p1_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
    ggplot2::geom_point() +
    ggpubr::stat_cor(label.y=lab*1.1)+
    ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
    facet_wrap(~variable, scales="free_x") +
    geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
    theme_classic()

p1_1

模块ZS值,网络属性,环境因子相关联合到一起

#网络模块和环境因子相关-网络指标和环境相关
dat.1 = env1 %>% rownames_to_column("id")
res = net.property.module.env (
  pst = pst,
  Top = 500,
  r.threshold= 0.8,
  p.threshold=0.05,
  env = dat.1,
  select.mod = c("model_1","model_2","model_3"),
  select.env = "env1")

#--模块和环境因子之间相关
p9 = res[[1]]
p9

#-网络属性和环境因子相关
p9 = res[[2]]
p9

时空组:网络相关分析

时空组:多组网络展示-网络性质稳定性等

这部分更新与2023年1月27日完成,可能会存在bug等,如有问题即使告知我解决;

这部分更新主要是为了解决一下几个问题: - 多个网络在一张图上需要调整映射颜色相同 - 分组并不是只有一列,很多时候有时间序列,有空间部位等,优化参数使得时间序列等多个分组的研究展示更加顺畅;

rm(list=ls())
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)
library(tidyfst)

ps.st = readRDS("./ps_TS.rds")
ps.st
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2432 taxa and 108 samples ]
## sample_data() Sample Data:       [ 108 samples by 4 sample variables ]
## tax_table()   Taxonomy Table:    [ 2432 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
## refseq()      DNAStringSet:      [ 2432 reference sequences ]
#时空组网络-分面网络图-解决填充颜色不一致问题#----

res = Facet.network (
    ps.st= ps.st,# phyloseq对象
    g1 = "Group",# 分组1
    g2 = "space",# 分组2
    g3 = "time",# 分组3
    ord.g1 = c("WT","KO","OE"),# 排序顺序
    ord.g2 = c("B","R") ,# 排序顺序
    ord.g3 = c("T1","T2","T3") ,# 排序顺序
    order = "time", # 出图每行代表的变量
    fill = "Phylum",
    size = "igraph.degree",
    layout_net = "model_maptree2",
    r.threshold=0.8,
    p.threshold=0.01,
    method = "spearman",
    select_layout = TRUE,
    clu_method = "cluster_fast_greedy",
    maxnode = 5
)
## [1] "B" "R"
## [1] "T1" "T2" "T3"
## [1] "WT" "KO" "OE"
p = res[[1]]

p

# ggsave("cs1.pdf",p,width = 5*5,height = 4*3)

#--时空组网络-网络性质#-------
dat = net.pro.ts(dat = res[[2]][[1]])
head(dat)
##                           B.T1.WT              B.T1.KO             
## num.edges(L)              "81"                 "96"                
## num.pos.edges             "41"                 "60"                
## num.neg.edges             "40"                 "36"                
## num.vertices(n)           "100"                "106"               
## Connectance(edge_density) "0.0163636363636364" "0.0172506738544474"
## average.degree(Average K) "1.62"               "1.81132075471698"  
##                           B.T1.OE              B.T2.WT             
## num.edges(L)              "66"                 "65"                
## num.pos.edges             "46"                 "35"                
## num.neg.edges             "20"                 "30"                
## num.vertices(n)           "96"                 "91"                
## Connectance(edge_density) "0.0144736842105263" "0.0158730158730159"
## average.degree(Average K) "1.375"              "1.42857142857143"  
##                           B.T2.KO              B.T2.OE             
## num.edges(L)              "106"                "78"                
## num.pos.edges             "67"                 "55"                
## num.neg.edges             "39"                 "23"                
## num.vertices(n)           "105"                "109"               
## Connectance(edge_density) "0.0194139194139194" "0.0132517838939857"
## average.degree(Average K) "2.01904761904762"   "1.43119266055046"  
##                           B.T3.WT              B.T3.KO             
## num.edges(L)              "73"                 "92"                
## num.pos.edges             "37"                 "55"                
## num.neg.edges             "36"                 "37"                
## num.vertices(n)           "94"                 "103"               
## Connectance(edge_density) "0.0167009837565774" "0.0175138016371597"
## average.degree(Average K) "1.5531914893617"    "1.78640776699029"  
##                           B.T3.OE              R.T1.WT             
## num.edges(L)              "77"                 "72"                
## num.pos.edges             "53"                 "36"                
## num.neg.edges             "24"                 "36"                
## num.vertices(n)           "105"                "96"                
## Connectance(edge_density) "0.0141025641025641" "0.0157894736842105"
## average.degree(Average K) "1.46666666666667"   "1.5"               
##                           R.T1.KO            R.T1.OE             
## num.edges(L)              "95"               "82"                
## num.pos.edges             "51"               "63"                
## num.neg.edges             "44"               "19"                
## num.vertices(n)           "107"              "108"               
## Connectance(edge_density) "0.01675189560924" "0.0141917618553133"
## average.degree(Average K) "1.77570093457944" "1.51851851851852"  
##                           R.T2.WT              R.T2.KO             
## num.edges(L)              "67"                 "107"               
## num.pos.edges             "33"                 "60"                
## num.neg.edges             "34"                 "47"                
## num.vertices(n)           "95"                 "106"               
## Connectance(edge_density) "0.0150055991041433" "0.0192273135669362"
## average.degree(Average K) "1.41052631578947"   "2.0188679245283"   
##                           R.T2.OE              R.T3.WT             
## num.edges(L)              "81"                 "72"                
## num.pos.edges             "60"                 "39"                
## num.neg.edges             "21"                 "33"                
## num.vertices(n)           "104"                "97"                
## Connectance(edge_density) "0.0151232262882748" "0.0154639175257732"
## average.degree(Average K) "1.55769230769231"   "1.48453608247423"  
##                           R.T3.KO              R.T3.OE             
## num.edges(L)              "94"                 "80"                
## num.pos.edges             "56"                 "57"                
## num.neg.edges             "38"                 "23"                
## num.vertices(n)           "99"                 "108"               
## Connectance(edge_density) "0.0193774479488765" "0.0138456213222568"
## average.degree(Average K) "1.8989898989899"    "1.48148148148148"
#--时空组的网络稳定性
#--时空组-稳定性1模块比对#--------
res = module.compare.m.ts(ps.st = ps.st, N = 200,
  degree = TRUE,zipi = FALSE,r.threshold= 0.8,
  p.threshold=0.05,method = "spearman",
  padj = F,n = 3,
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  zoom = 0.3,# 控制小圈大小
  b.x = 1.3,
  b.y = 1.3)
## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 2 8
## [1] 2 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 2 8
## [1] 2 8
## NULL
## [1] 1 8
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## [1] 1 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 1 8
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 1 8
## [1] 1 8
## [1] 2 8
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
p = res[[1]]
p

#提取数据
dat = res[[2]]
#-作图数据提取
dat = res[[3]]

#时空组-稳定性2 鲁棒性随即取出#-----
res = Robustness.Random.removal.ts(
  ps.st = ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",
  g2 = "space",
  g3 = "time")

#
res[[1]][[1]]

res[[1]][[2]]

#时空组-稳定性2-2 鲁棒性随即取出#-----
res = Robustness.Targeted.removal.ts(
  ps.st,
  N = 200,
  degree = TRUE,
  zipi = FALSE,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time")

res[[1]][[1]]

res[[1]][[2]]

#时空组-稳定性3负相关比例#-------
res = negative.correlation.ratio.ts(
  ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 =NULL,# 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL# 排序顺序
)
# 导出图片
res[[1]]

# 导出数据
dat = res[[2]]
head(dat)
##        ID    ratio group time space label   Group
## 1 R.T1.KO 46.31579    KO   T1     R  T1.R R.T1.KO
## 2 R.T1.OE 23.17073    OE   T1     R  T1.R R.T1.OE
## 3 R.T1.WT 50.00000    WT   T1     R  T1.R R.T1.WT
## 4 R.T2.KO 43.92523    KO   T2     R  T2.R R.T2.KO
## 5 R.T2.OE 25.92593    OE   T2     R  T2.R R.T2.OE
## 6 R.T2.WT 50.74627    WT   T2     R  T2.R R.T2.WT
#时空组-稳定性4 网络易损性#------
library("pulsar")
res = natural.con.microp.ts(
  ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 =NULL,# 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL,# 排序顺序
  norm = F,
  end = N -2,
  start = 0,
  con.method = "pulsar"

)

res[[1]]

#时空组-稳定性5 组成稳定性#-------
res = community.stability.ts(
  ps.st = ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "space",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  map.art = NULL, # 人工输入的分组 默认为NULL
  time = F,# 稳定性是否有时间序列
  ord.map = TRUE# map文件是否是已经按照pair要求进行了排序
)

res[[1]]

# res[[2]]

时空组:多网络单独绘制却填充相同颜色

#---时空组双网络手动填充相同颜色#-------
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)

ps.st = readRDS("./ps_TS.rds")
ps.st
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2432 taxa and 108 samples ]
## sample_data() Sample Data:       [ 108 samples by 4 sample variables ]
## tax_table()   Taxonomy Table:    [ 2432 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
## refseq()      DNAStringSet:      [ 2432 reference sequences ]
#-----第一组网络绘制

res = Facet.network(
  ps.st= ps.st,# phyloseq对象
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 = c("WT","KO","OE"),# 排序顺序
  ord.g2 = c("B","R") ,# 排序顺序
  ord.g3 = c("T1","T2","T3") ,# 排序顺序
  order = "time", # 出图每行代表的变量
  fill = "Phylum",
  size = "igraph.degree",
  layout_net = "model_maptree2",
  r.threshold=0.8,
  p.threshold=0.01,
  method = "spearman",
  select_layout = TRUE,
  clu_method = "cluster_fast_greedy",
  maxnode = 5
)
## [1] "B" "R"
## [1] "T1" "T2" "T3"
## [1] "WT" "KO" "OE"
p = res[[1]]
p

#--指定颜色映射,为了多个图使用同一个颜色#-------



#--设定的参数一致
fill = "Phylum"
size = "igraph.degree"
maxnode = 5
row.num = 6

#-从结果中提取网络节点和边数据
node  = res[[2]][[2]]
head(node)
##         ID         X1        X2 elements igraph.degree igraph.closeness
## 1  ASV_101  1.8922077 -16.17873  ASV_101             3                1
## 2 ASV_1050 10.4488008  11.60225 ASV_1050             1                1
## 3  ASV_113 -0.3111552  18.45045  ASV_113             0                0
## 4 ASV_1148 11.7082072 -13.80247 ASV_1148             0                0
## 5 ASV_1149  2.1328120  18.18372 ASV_1149             1                1
## 6  ASV_115  7.1706596 -16.75757  ASV_115             0                0
##   igraph.betweenness igraph.cen.degree   group  Kingdom         Phylum
## 1                  0                 3 B.T1.WT Bacteria Proteobacteria
## 2                  0                 1 B.T1.WT Bacteria Proteobacteria
## 3                  0                 0 B.T1.WT Bacteria Actinobacteria
## 4                  0                 0 B.T1.WT Bacteria     Firmicutes
## 5                  0                 1 B.T1.WT Bacteria     Firmicutes
## 6                  0                 0 B.T1.WT Bacteria Proteobacteria
##                 Class            Order             Family          Genus
## 1          Unassigned       Unassigned         Unassigned     Unassigned
## 2 Alphaproteobacteria Sphingomonadales  Sphingomonadaceae   Sphingomonas
## 3      Actinobacteria  Actinomycetales  Microbacteriaceae      Agromyces
## 4             Bacilli       Bacillales     Planococcaceae Lysinibacillus
## 5             Bacilli       Bacillales Paenibacillaceae_1     Unassigned
## 6  Betaproteobacteria  Burkholderiales   Oxalobacteraceae     Unassigned
##                      Species Group time space links nodes
## 1                 Unassigned    WT   T1     B   162   100
## 2 Sphingomonas_daechungensis    WT   T1     B   162   100
## 3                 Unassigned    WT   T1     B   162   100
## 4                 Unassigned    WT   T1     B   162   100
## 5                 Unassigned    WT   T1     B   162   100
## 6                 Unassigned    WT   T1     B   162   100
##                               label
## 1 B.T1.WT: (nodes: 100; links: 162)
## 2 B.T1.WT: (nodes: 100; links: 162)
## 3 B.T1.WT: (nodes: 100; links: 162)
## 4 B.T1.WT: (nodes: 100; links: 162)
## 5 B.T1.WT: (nodes: 100; links: 162)
## 6 B.T1.WT: (nodes: 100; links: 162)
edge = res[[2]][[3]]
head(edge)
##         X2       Y2  OTU_2   OTU_1 weight       X1       Y1 cor   group Group
## 1 6.251566 6.418839 ASV_30  ASV_21     -1 5.453796 4.025147   - R.T3.OE    OE
## 2 6.251566 6.418839 ASV_30  ASV_29      1 4.842470 7.579802   + R.T3.OE    OE
## 3 6.251566 6.418839 ASV_30  ASV_55      1 3.675405 8.282087   + R.T3.OE    OE
## 4 6.251566 6.418839 ASV_30 ASV_290      1 3.529926 5.996121   + R.T3.OE    OE
## 5 4.842470 7.579802 ASV_29  ASV_21     -1 5.453796 4.025147   - R.T3.OE    OE
## 6 4.842470 7.579802 ASV_29  ASV_30      1 6.251566 6.418839   + R.T3.OE    OE
##   time space links nodes                             label
## 1   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
## 2   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
## 3   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
## 4   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
## 5   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
## 6   T3     R   160   108 R.T3.OE: (nodes: 108; links: 160)
cb_palette <- c("#ed1299", "#09f9f5", "#246b93", "#cc8e12", "#d561dd", "#c93f00", "#ddd53e",
                "#4aef7b", "#e86502", "#9ed84e", "#39ba30", "#6ad157", "#8249aa", "#99db27", "#e07233", "#ff523f",
                "#ce2523", "#f7aa5d", "#cebb10", "#03827f", "#931635", "#373bbf", "#a1ce4c", "#ef3bb6", "#d66551",
                "#1a918f", "#ff66fc", "#2927c4", "#7149af" ,"#57e559" ,"#8e3af4" ,"#f9a270" ,"#22547f", "#db5e92",
                "#edd05e", "#6f25e8", "#0dbc21", "#280f7a", "#6373ed", "#5b910f" ,"#7b34c1" ,"#0cf29a" ,"#d80fc1",
                "#dd27ce", "#07a301", "#167275", "#391c82", "#2baeb5","#925bea", "#63ff4f")
tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf = data.frame(id = tem1,color =cb_palette[1:length(tem1)] )
head(tabf)
##                id   color
## 1  Proteobacteria #ed1299
## 2  Actinobacteria #09f9f5
## 3      Firmicutes #246b93
## 4   Bacteroidetes #cc8e12
## 5      Unassigned #d561dd
## 6 Verrucomicrobia #c93f00
colnames(tabf)[1] = fill
node1 = node %>% left_join(tabf,by = fill)

p2 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
                                  color = cor),
                              data = edge, size = 0.03,alpha = 0.5) +
  geom_point(aes(X1, X2,
                 size = !!sym(size),
                 fill = color ),
             pch = 21, data = node1,color = "gray40") +
  facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  # geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  scale_colour_manual(values = c("#6D98B5","#D48852")) +
  scale_fill_manual(values  = tabf$color,labels = tabf[,fill]) +
  scale_size(range = c(0.8, maxnode)) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)
  ) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()
  ) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2

# ggsave("cs1.pdf",p2,width = 5*5,height = 4*3)



#-------另外一组网络绘制



data(ps)

res2 = Facet.network (
  ps.st= ps,# phyloseq对象
  g1 = "Group",# 分组1
  # g2 = "space",# 分组2
  # g3 = "time",# 分组3
  ord.g1 = c("WT","KO","OE"),# 排序顺序
  # ord.g2 = c("B","R") ,# 排序顺序
  # ord.g3 = c("T1","T2","T3") ,# 排序顺序
  order = "time", # 出图每行代表的变量
  fill = "Phylum",
  size = "igraph.degree",
  layout_net = "model_maptree2",
  r.threshold=0.8,
  p.threshold=0.01,

  method = "spearman",
  select_layout = TRUE,
  clu_method = "cluster_fast_greedy",
  maxnode = 5
)
## [1] "WT" "KO" "OE"
p3 = res2[[1]]
p3

#-提取节点和边数据
node  = res2[[2]][[2]]
head(node)
##        ID        X1           X2 elements igraph.degree igraph.closeness
## 1   ASV_1  6.394069  -0.02808431    ASV_1             1                1
## 2  ASV_10  3.652617 -21.02187428   ASV_10             0                0
## 3 ASV_100  7.078794  10.35742308  ASV_100             1                1
## 4 ASV_101 -3.792132  16.03963922  ASV_101             1                1
## 5 ASV_102  8.292658  -9.25473746  ASV_102             1                1
## 6 ASV_103  4.025401  21.04072916  ASV_103             0                0
##   igraph.betweenness igraph.cen.degree group  Kingdom         Phylum
## 1                  0                 1  ..WT Bacteria Actinobacteria
## 2                  0                 0  ..WT Bacteria Proteobacteria
## 3                  0                 1  ..WT Bacteria Proteobacteria
## 4                  0                 1  ..WT Bacteria Proteobacteria
## 5                  0                 1  ..WT Bacteria Proteobacteria
## 6                  0                 0  ..WT Bacteria Proteobacteria
##                 Class           Order              Family           Genus
## 1      Actinobacteria Actinomycetales Thermomonosporaceae      Unassigned
## 2 Alphaproteobacteria     Rhizobiales        Rhizobiaceae       Rhizobium
## 3  Betaproteobacteria Burkholderiales      Comamonadaceae  Hydrogenophaga
## 4          Unassigned      Unassigned          Unassigned      Unassigned
## 5 Gammaproteobacteria      Unassigned          Unassigned      Unassigned
## 6 Alphaproteobacteria     Rhizobiales  Phyllobacteriaceae Phyllobacterium
##                        Species Group time space links nodes
## 1                   Unassigned    WT              138   105
## 2                   Unassigned    WT              138   105
## 3    Hydrogenophaga_intermedia    WT              138   105
## 4                   Unassigned    WT              138   105
## 5                   Unassigned    WT              138   105
## 6 Phyllobacterium_bourgognense    WT              138   105
##                            label
## 1 ..WT: (nodes: 105; links: 138)
## 2 ..WT: (nodes: 105; links: 138)
## 3 ..WT: (nodes: 105; links: 138)
## 4 ..WT: (nodes: 105; links: 138)
## 5 ..WT: (nodes: 105; links: 138)
## 6 ..WT: (nodes: 105; links: 138)
edge = res2[[2]][[3]]
head(edge)
##           X2         Y2  OTU_2   OTU_1 weight         X1         Y1 cor group
## 1  1.9872552  1.5012669 ASV_22  ASV_32      1 -3.0932354 -2.3169033   +  ..OE
## 2  1.9872552  1.5012669 ASV_22  ASV_77      1 -0.4037329 -0.8876355   +  ..OE
## 3  1.9872552  1.5012669 ASV_22  ASV_95     -1 -0.9859822  2.3370501   -  ..OE
## 4  1.9872552  1.5012669 ASV_22 ASV_214      1  2.3530761 -1.0455685   +  ..OE
## 5  1.9872552  1.5012669 ASV_22 ASV_156      1 -3.4791850  0.5142463   +  ..OE
## 6 -0.4037329 -0.8876355 ASV_77  ASV_32      1 -3.0932354 -2.3169033   +  ..OE
##   Group time space links nodes                         label
## 1    OE              148    96 ..OE: (nodes: 96; links: 148)
## 2    OE              148    96 ..OE: (nodes: 96; links: 148)
## 3    OE              148    96 ..OE: (nodes: 96; links: 148)
## 4    OE              148    96 ..OE: (nodes: 96; links: 148)
## 5    OE              148    96 ..OE: (nodes: 96; links: 148)
## 6    OE              148    96 ..OE: (nodes: 96; links: 148)
node$group %>% unique()
## [1] "..WT" "..KO" "..OE"
#----第二组网络用第一组颜色填充#---------


tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf2 = data.frame(id = tem1 )
colnames(tabf2)[1] = fill
head(tabf2)
##           Phylum
## 1 Actinobacteria
## 2 Proteobacteria
## 3  Bacteroidetes
## 4     Unassigned
## 5    Chloroflexi
## 6     Firmicutes
# tabf2[1,1] = "wentao"
tabf3 = tabf2 %>% left_join(tabf,by = fill)

num = tabf3$color[is.na(tabf3$color)] %>% length()
if (num == 0) {
  tabf.f = tabf3
}else{
  tem = tabf3 %>% filter(is.na(color))
  tem$color = setdiff(cb_palette,tabf$color)[1:length(tem$color)]
  tabf.f = rbind(tabf3,tem)
}

colnames(tabf.f)[1] = fill
node1 = node %>% left_join(tabf.f,by = fill)

head(node1)
##        ID        X1           X2 elements igraph.degree igraph.closeness
## 1   ASV_1  6.394069  -0.02808431    ASV_1             1                1
## 2  ASV_10  3.652617 -21.02187428   ASV_10             0                0
## 3 ASV_100  7.078794  10.35742308  ASV_100             1                1
## 4 ASV_101 -3.792132  16.03963922  ASV_101             1                1
## 5 ASV_102  8.292658  -9.25473746  ASV_102             1                1
## 6 ASV_103  4.025401  21.04072916  ASV_103             0                0
##   igraph.betweenness igraph.cen.degree group  Kingdom         Phylum
## 1                  0                 1  ..WT Bacteria Actinobacteria
## 2                  0                 0  ..WT Bacteria Proteobacteria
## 3                  0                 1  ..WT Bacteria Proteobacteria
## 4                  0                 1  ..WT Bacteria Proteobacteria
## 5                  0                 1  ..WT Bacteria Proteobacteria
## 6                  0                 0  ..WT Bacteria Proteobacteria
##                 Class           Order              Family           Genus
## 1      Actinobacteria Actinomycetales Thermomonosporaceae      Unassigned
## 2 Alphaproteobacteria     Rhizobiales        Rhizobiaceae       Rhizobium
## 3  Betaproteobacteria Burkholderiales      Comamonadaceae  Hydrogenophaga
## 4          Unassigned      Unassigned          Unassigned      Unassigned
## 5 Gammaproteobacteria      Unassigned          Unassigned      Unassigned
## 6 Alphaproteobacteria     Rhizobiales  Phyllobacteriaceae Phyllobacterium
##                        Species Group time space links nodes
## 1                   Unassigned    WT              138   105
## 2                   Unassigned    WT              138   105
## 3    Hydrogenophaga_intermedia    WT              138   105
## 4                   Unassigned    WT              138   105
## 5                   Unassigned    WT              138   105
## 6 Phyllobacterium_bourgognense    WT              138   105
##                            label   color
## 1 ..WT: (nodes: 105; links: 138) #09f9f5
## 2 ..WT: (nodes: 105; links: 138) #ed1299
## 3 ..WT: (nodes: 105; links: 138) #ed1299
## 4 ..WT: (nodes: 105; links: 138) #ed1299
## 5 ..WT: (nodes: 105; links: 138) #ed1299
## 6 ..WT: (nodes: 105; links: 138) #ed1299
p4 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
                                  color = cor),
                              data = edge, size = 0.03,alpha = 0.5) +
  geom_point(aes(X1, X2,
                 fill =color,
                 size = !!sym(size)),
             pch = 21, data = node1,color = "gray40") +
  facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  # geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  scale_colour_manual(values = c("#6D98B5","#D48852")) +
  scale_fill_manual(values  = tabf.f$color,labels = tabf.f[,fill]) +
  scale_size(range = c(0.8, maxnode)) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)
  ) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()
  ) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
# guides(fill=FALSE)
p4

# ggsave("cs2.pdf",p4,width = 5*3,height = 4*1)
library(patchwork)
p2|p4